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We formulate a quantum waveguide theory of the Josephson effect in multiband superconductors, 
with special emphasis on iron-based materials. By generalizing the boundary conditions of the 
scattering problem, we first determine the Andreev levels spectrum and then derive an explicit 
expression for the Josephson current which generalizes the formula of the single band case. In 
deriving the results, we provide a second quantization field theory allowing to evaluate the current- 
phase relation and the Josephson current fluctuations in multiband systems. We present results 
for two different order parameter symmetries, namely s± and s++, which are relevant in multiband 
systems. The obtained results show that the s± symmetry can support n states which are absent in 
the s++ case. We also argue a certain fragility of the Josephson current against phase fluctuations in 
the s++ case. The temperature dependence of the Josephson critical current is also analyzed and we 
find, for both the order parameter symmetries, remarkable violations of the Ambegaokar-Baratoff 
relation. The results are relevant in view of possible experiments aimed at investigating the order 
parameter symmetry of multiband superconductors using mesoscopic Josephson junctions. 

PACS numbers: 74.45.+C, 74.50.+r, 74.70.Xa, 85.25.Cp 


I. INTRODUCTION 

Multiband superconductivity [I] has been theoretically 
suggested few years later the BCS formulation of the su¬ 
perconducting state [1|. The discovery of the supercon¬ 
ductivity in MgB 2 0 in 2001 and the observation of su¬ 
perconductivity in Fe-based compounds (FeBS) SI in 
2008, renewed the interest towards multiband supercon¬ 
ductivity. One main reason is the comprehension of the 
pairing symmetry in this class of superconductors which 
is still object of intensive investigation. In the case of 
MgB 2 there is a clear consensus towards the picture of 
two coexisting in-phase superconducting gaps (s++ pair¬ 
ing) jg], while in FeBS the experimental evidences seem 
to be favorable, in some cases, to the s± pairing, implying 
that both the electron-like and the hole-like band develop 
an s-wave superconducting state with order parameters 
of opposite sign [7]. Thus, in the symmetry, the super¬ 
conducting gap exhibits a sign reversal between a and /3 

bands which is absent in the s_|_|_ symmetry. This phase 

difference can be probed using point contact Andreev 
reflection spectroscopy (PCARS) which is considered as 
one of the high-resolution phase-sensitive techniques to 
investigate the superconducting order parameter [8Hll|. 
This technique has recently been applied to gain insight 
on the properties of the FeBS [III El- However, in con¬ 
trast to the high-T), cuprate superconductivity manifest¬ 
ing a d-wave symmetry, phase sensitive experiments ap- 
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pear difficult in FeBS since the s++ and s± pairings have 
the same crystallographic symmetry and, as a matter of 
fact, so far no experiment has been decisive in discrim¬ 
inating between the two symmetries. A complementary 
tool to gain information on the symmetry of the order 
parameter is the current-phase relation of the Josephson 
current. In fact, along with PCARS, the Josephson effect 
has been used in the past as a probe of electronic prop¬ 
erties in superconductors, included the order parameter 

symmetry (e.g. (M3)- 

In this paper, in order to discriminate between differ¬ 
ent symmetries of the order parameter (s++ and s±), we 
develop a quantum waveguides approach for the Joseph¬ 
son effect in a multiband superconducting junction. The 
presence of more than one band in the superconduc¬ 
tor implies that extra scattering channels are present 
at the interface, a physical situation which is analo¬ 
gous to a quantum waveguides theory problem mm 
as recently suggested by PCARS modeling of a normal- 
multiband superconductor junction [j20| • We evaluate the 
Josephson current carried by Andreev bound levels and 
demonstrate that several distinctive features of the s± 
and s ++ symmetries of the order parameter can be ev¬ 
idenced in the current-phase relation. The results are 
particular relevant with respect to the fast progresses 
in nanofabrication techniques plj which allow now to 
explore the Josephson effect in the mesoscopic junction 
regime, where supercurrent flows through a small num¬ 
ber of channels. In this kind of nanometric junctions the 
effect of a different symmetry, s± or s++, should emerge 
and the experimental results easily compared with the 
basic theory of the Josephson effect in these materials. 
The paper is organized as follows. In Section [TT] we 
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present the quantum waveguide model of a multiband 
superconductor coupled to another identical multiband 
superconductor (symmetrical junction) and derive the 
spectral equation of the Andreev bound states. In Sec¬ 
tion m we calculate the Josephson current by using a 
field theory formalism in second quantization. In Section 
m we discuss the results for the current-phase relation 
and the temperature dependence of the critical current. 
In Section 0 we draw the conclusions. Details on the 
computation of the scattering coefficients are reported 
in Appendix [A] The theory of the magnetic field depen¬ 
dence of the critical current of a Josephson junction is 
briefly recalled in Appendix iBl 



FIG. 1. Quantum waveguide schematization of a two-band 
Josephson junction. Each branch represents a superconduct¬ 
ing band, a or /?> in the left ( Sl ) or right (Sr) electrode. 
The wave functions of supeconducting bands are given by 
'& S a L p R ( x )- The node at x = 0 represents the interface be¬ 
tween the left and right superconducting banks. 


II. MODEL AND THEORY 

We formulate a ballistic theory of the Josephson effect 
describing multiband superconducting junctions. The 
theory allows to consider an arbitrary number of bands, 
which are treated as network branches of an effective 
quantum waveguide model. The proposed approach 
could be applied to the MgB 2 as well as the FeBS case. 
In order to develop the theory, hereafter we refer to the 
specific case of FeBS for which the s ++ and s± sym¬ 
metries have been suggested. The Josephson effect in 
FeBS has been already studied by using a number of dif¬ 
ferent methods [22 h 41]1. However, most of these stud¬ 
ies, except few cases [25l-l27j]. deal with the so-called hy¬ 
brid Josephson junction in which the junction is formed 
by a conventional (s-wave) superconductor and a (s± or 
s++) FeBS. We specialize here to the case of an all- 
FeBS coplanar Josephson junction, in which both the 
electrodes are FeBS materials (symmetric junction), as 
occurs for instance in grain boundary junctions. So far, 
several Josephson junctions using thin films have been 
fabricated |12ho l| on bicrystal substrates (for a review 
on Fe-based Josephson junctions, see Ref.|52j] and ref¬ 
erences therein) and thus a theoretical effort along this 
direction is needed. 

Thus we consider a Josephson junction model in which 


the electrodes are two-band FeBS. To this end the junc¬ 
tion is represented by a network of four one-dimensional 
branches connected to a single node point x = 0 (see 
Fig. m, x being the coordinate along the propagation 
direction normal to the interface. Each superconduct¬ 
ing branch represents the i-th band on the left or the 
right side of the junction, while the band wave func¬ 
tion Hq = 9(—x)^ L (x) + 9(x)' i ^f R (x), in the clean limit, 
obeys the equation 

Hj(x) A j(x) 

A*(x) -H*(x) 

where j € { a , /3} is the band index, while 

Aj (x) = Aje^Oi-x) + A^e^+^eix), (2) 

(9(x) = 0, for x < 0, 9(x) = 1 for x > 0) are the two 
coexisting pair potentials. The operator Hj(x ) represents 
the single particle Bogoliubov-de Gennes Hamiltonian in 
the two bands, which reads 


'I', = (1) 


H 0 (x) = 


h 2 d 2 
2rrij dx 2 


E f . 


( 3 ) 


In writing Eq. m we are neglecting the proximity effect 
and assuming that the inhomogeneous character of the 
gap in the junction is captured by choosing A a = A^ = 0 
just at the node (short junction). The two gap values 
A a , Ap, with A q < Ap, are assumed to be the same 
in the two superconductive leads. The quantity ip is the 
gauge invariant phase difference between the two super¬ 
conductive regions, ip a and ipp the internal pair potential 
phases. In the case of s±-wave gap model, ipp — ip a = 7r 
and the two gap have opposite sign, while in the stan¬ 
dard two band model, with same sign gaps, ipp — ip a = 0. 
The quantities nij in Eq. m are the effective masses of 
quasiparticles in the j-th superconducting branch and 
are material-dependent quantities. We also introduce 
a single-particle node potential U(x) which is different 
from zero only for x = 0 and can be modeled as the 
usual Blonder-Tinkham-Klapwijk [53] interface potential 
U(x) = UoS(x), even though this limiting assumption 
is not required to develop the theory. The node poten¬ 
tial U(x) allows the modeling of a FeBS/insulator/FeBS 
or FeBS/normal-metal/FeBS junction in which a nor¬ 
mal scattering at the interface reduces the transparency 
of the junction. The junction barrier strength can be 
still characterized by introducing a Blonder-Tinkham- 
Klapwijk dimensionless parameter Z = mUo/(h 2 k F ) (jn 
being the bare electron mass and kp = 2 mE F /h 2 ) within 
the boundary conditions of the scattering problem. The 
modified boundary conditions used here (see Eqs.©, ©) 
account for band-sensitive scattering effects and have 
been already introduced in Ref. [20( to describe the dif¬ 
ferential conductance of a normal-metal/FeBS junction; 
in the following we discuss their generalization to the 
Josephson junction case. The potentials U(x) and A j(x) 
are responsible for the normal scattering and the scat¬ 
tering of electrons into holes (Andreev scattering) at the 
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interface, respectively. The four wave functions, 'I'J 7 ' (x), 
'I'| R (*) (j = a, ft), one for each branch in the two elec¬ 
trodes, can be written in terms of the eigenstates of the 
local Hamiltonians as 


4> j{x ) = 


uf (x) 
vhx) 


a, 


Uje ZV!j 


e iPhX 


Vje llpi 


(4) 


for x < 0 or 


*?*(*) = 


u^( x) 
Vj(x) 


Cn 


u 

v^e 


3 ,1 rl . 

,^ -^iVj+v) ' e + «7 


3 ■/ , 1 e~ lp h x 

Uje -<vi+<p) I e 


(5) 


for x > 0. Here u a and v a , and up and vp are the 
Bogoliubov coefficients for the first and second band, re¬ 
spectively 



'1 

A n,-\ 1 

1/2 

ri 

/ HAl 

Uj = 

2 


,Vj = 

[2 



( 6 ) 


with Qj = Aj(T) 2 — E 2 , j G { a , ft}, while T represents 
the temperature of the thermal bath. The wave vectors 

I 


of the two bands 

Pi,h = \Jr]k 2 F ±2im^ J /h 2 
are well approximated by the expressions: 

Pe,h - r a k F (l ± , (7) 

Pe,h - r P k F ± ix\! A/3 2 /A« 2 - E 2 /A a 2 ^j , 

which are valid under the assumptions Ep A^, while 
the different signs ” refer to electronlike or hole¬ 

like excitations respectively. The coefficient r 2 = rrij/m 
represents the effective mass of the j-th band mea¬ 
sured in units of the bare electron mass to. In writ¬ 
ing Eqs. © we have introduced the dimensionless factor 
X = A a (T)/2E F or x = A a (T)/(k F v F h ), which repre¬ 
sents the ratio between the Fermi wavelength k F 1 and the 
coherence length £(T) ~ hv F /A a (T). Its zero tempera¬ 
ture value, xo, can be also defined as xo = A a (0)/(2E F ). 
Since we are interested in localized subgap states, in Eq. 
( 13 ) we assume E < A Q (0), ensuring that the wave func¬ 
tions 4 , J I ‘(:r), 4/^ R ( x) (j = a, ft) decay exponentially for 
|a:| —>• oo. 

The quasiparticle wave functions obey the generalized 
matching quantum waveguide conditions |20l 




Tf(0) = ^4-^(0), *£*(0) = stf^(0) 


( 8 ) 


d_ 

dx 


'a 'ft 


x=0 


d_ 

dx 


Equations © and © guarantee the conservation of the 
charge current at the node * = 0 (quantum Kirchhoff’s 
law) and generalize the waveguide boundary conditions 
given in Ref. |2(j to the Josephson junction case. The 
three new parameters s a ,sp,s characterize the inter¬ 
face and make the scattering and the current flow band- 
sensitive. Indeed, the overlap between the wave functions 
on different sides of the junction may favour scattering 
events towards a specific band, a or ft. The existence 
of these ’band coupling parameters’ implies that a dis¬ 
continuity of the wave functions may occur at the node 
(x = 0) when at least one of the coupling parameters 
SonSp,s is different from one. The meaning of this dis¬ 
continuity has been discussed in Ref. j2(j in connection 
with the theory of PCARS in multiband superconduc¬ 
tors. The physical meaning of the overlap parameters 
s a ,sp,s can be easily understood observing that setting, 
for instance, s a = 0 implies 4^(0) = 0. Thus the wave 
function presents a vanishing overlap at the in¬ 

terface with the remaining branch wave functions. As a 


Is 1 

— vl/ 0 ^ _i_ s — \jy 

r 2 a w * 2 * 

' cv ' R 


S L 


- x=0 


2k F Z^ L (0). 


(9) 


consequence, in this case, the band a on the right side of 
the junction is excluded from the transport. The latter 
point can be further clarified by looking at the prefac¬ 
tor of the first term of Eq. ©, namely s a (m a /m) _1 , 
which governs the particle flow through the a band on 
the right side of the junction. Lowering s a has the same 
effect of increasing m a , causing a lowering of the particles 
flux through the considered band. In the equal coupling 
case (s Q = S /3 = s = l) the partitioning of the current 
among different network branches is entirely governed by 
the bulk mass ratios Vj and only weakly affected by the 
interface scattering potential strength Z. Thus the pa¬ 
rameters s a ,sp,s can be properly interpreted as overlap 
or coupling factors which define the fraction of current 
flowing through a specific band. 

Eqs. ©-© provide a linear system of equations for the 
eight unknown coefficients a,j, bj , Cj, dj ( j € {a, ft}), 
appearing in the wave functions (see Appendix [Aj. The 
bound state spectrum can be determined imposing the 
condition that this homogeneous system of equations has 
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a nontrivial solution. Within the Andreev’s approxima- critical condition is written as 
tion, i.e. p“ ~ p^ ~ r a kp and p@ ~ p^ ~ rpkp, this 


S |(e 2 - 7 2 )(1 - e 2 - 2s 2 e 2 - 4(e 2 - 1) + 2 S 2 cos*?) 

r 2 
1 a 

s 2 (e 2 — l)(2ss^e 2 + (s 2 + s^)(e 2 — 7 2 ) — 2ss 2 /j 'y 2 cost/?) 

12 h 

r P 

2sspVe 2 - ly/e 2 - 7 2 (s((l + s 2 )e 2 - e I,5 7) + s|((l + s 2 )e 2 - e l5 s 2 7) - e l<5 (ss 2 + sfy'ycospj 


4 S 2 Z 2 (e 2 -l)( e 2 - 7 2 )=0, 


( 10 ) 


where 6 = ipp — ip a , 7 = Ap/A a and e = E/A a (T). 
Choosing S = 0 or 6 = 7 r reflects the internal pairing sym¬ 
metry s ++ or s-t, respectively. Two sets of energy levels 
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FIG. 2. Bound states energy levels E + (ip) and E~(ip) (nor¬ 
malized to A a(T)) calculated through Eq. (1101) . Panels (a), 
(c): s± symmetry. Panels (b), (d): s++ symmetry. Z and 7 
have been fixed to Z = 0.03 and 7 = 3, the various curves are 
obtained by letting parameter s a to vary in the range (0.3,1) 
in steps of 0.1, while sp = s = 1. Panels (a) and (b): r a = 1.8, 
rp = 1. Panels (c) and (d): r a = 1, rp = 1.8, the blue line 
curve in panel (c) is the Andreev level of a conventional s-wave 
junction with rjZ = 0.03. 


E + (ip) and E (ip) are found in solving Eq. m, which 
correspond to electronlike and holelike quasiparticles, re¬ 
spectively. One finds that the Andreev bound states 
never appear beyond A a . When s a = sp = s = 1, in 
the two limits r a —> 00 or rp —> 00 (exclusion of the a or 
/3 band respectively) Eq. (ITOl) provides the well known s- 


wave result E = Eb = Aj (cos 2 (ip/ 2) + Z 2 )/(l + Z 2 ), 

with Zj = rjZ {54j , j being a or /3 (see Fig. [2] (c)). 

In Fig. [5] we show Andreev bound states spectra com¬ 
puted using Eq. m- Once the mass ratios r :j have been 
fixed, different Josephson couplings can be obtained de¬ 
pending on the kind of the order parameter symmetry 
and on the choice of the overlap factors. To better ex¬ 
plain this point, we present different curves obtained by 
letting the parameter s a to vary in the range (0.3,1) in 
steps of 0.1 while sp = s = 1. This choice of the overlap 
factors corresponds to weaken progressively the coupling 
of the band \I/f K from the remaining bands. Panels (a) 
and (b) ( r a = 1.8, rp = 1) show a generalized lowering of 
the levels toward the zero energy state in the case of s± 
symmetry (panel (a)), which is absent in the case of s_|_+ 
symmetry (panel (b)) and represents the main difference 
between the two symmetries in this case. A further differ¬ 
ence is illustrated in panels (c) and (d) ( r a = 1, rp = 1.8). 
In the case of s± symmetry (panel (c)) the energy lev¬ 
els undergo a gradual concavity change as the parameter 
s a is decreased. This concavity change is responsible for 
a 7r-shift in the current-phase relation, which is not ob¬ 
served in the s ++ case (Fig.[2])d)). Thus, as also observed 
elsewhere (25[, a 0 — n transition may be a distinctive 
signature of the s± symmetry. This point will be further 
explained in Section m where the current-phase rela¬ 
tions and the temperature dependence of the maximum 
Josephson current will be presented. 


III. JOSEPHSON CURRENT 


The Josephson current flowing through the junction 
could be directly computed from the Andreev bound 
states according to the formula [55| : 


Ij(p) « 


H£dAM tanh [AMl. 

H dip L 2 ksT - 


( 11 ) 


Due to the special character of the bound state of a multi¬ 
band junction, in the following we will characterize the 
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Josephson current by using a field theory formalism in 
second quantization. We verified that the latter proce¬ 
dure provides results consistent with those obtained using 

Eq. (fHD . 

The Cooper pair charge flow J c h through the junction 
can be computed as quantum-statistical average of the 
current density operator according to the expression 
Jch = where jj = ! ^(3x-^dx) rep¬ 

resents the first-quantization band operator [50]. The 
second quantized fields ipj a (V’j cr ) represent the annihi¬ 
lation (creation) operators of a spin a electron in the 
j -th band and obey a fermionic algebra. These fields 
provide an appropriate basis in the absence of supercon¬ 
ducting correlations. However, when the superconduc¬ 
tivity is established, i/jja are not eigenfields of the prob¬ 
lem and, for this reason, the Nambu representation can 
be conveniently used ^(z) = Eje{a,/3} |j> ® (V’jf ,$4.)*- 
The Nambu field (rc) is a non-local quantity describ¬ 
ing particle-hole excitations in different network branches 
and can be expanded in eigenfields of the Hamiltonian 
problem on the network. Limiting the expansion to the 
eigenfields 7^ describing the low-energy (sub-gap) states 
with energy ±E B G [—A q , A a ], we get: 

*(®) « 7 te~ iEBt/ ^ B (x) + 7 \e iE ° t/h ^ B {x), (12) 


where ip B ( x ) = 6(—x)ip^\x) + 9{x)ip^\x) is the wave 
function of the electron-like bound state having en¬ 
ergy eigenvalue E B > 0, while ^ B {x) = [Ej |j)(j| <8> 
ia y C\il> B (x) represents its time-reversed state associated 
with a hole-like state with energy —E B . The electron¬ 
like bound state is localized at the interface and extends 
over all the waveguide branches. Thus eigenstates of the 
local branch Hamiltonians can be used to expand ip B (x ), 
using the following decomposition: 


B \x)= y IJ>® 

ie{«,/3} 


Uj (x) ' 
Vj( x) J ’ 


(13) 


with v € [L, R] an index identifying the left (x < 0) or 
right (x > 0) side of the junction. Once ^(a:) has been 
expanded in eigenfields j a it is possible to recognize the 
fermionic fields ipjf in the expression: 


i-t « 7t e- iEBt/h Uj(x) 7+ e iE - t/n v^r, (14) 


the current conservation we have the freedom to evaluate 
the current in x = 0 + , and thus starting from Eq. (fl5l) 
we get (e = —|e|) (58j : 


J ch = —2evp Y - Lj ~ L (l c j| 2 - Mil 2 )tanh 

T j 

J J 

where the coefficients Cj and dj are calculated at the 
energy E B = e B (ip)A a (T), while e B (ip) is solution of 
Eq. (HUD- It’s worth mentioning here that the Fermi ve¬ 
locity in the «-th band is given by v^) = vf/Ti , while 
the quantity vf = Hkp/m is just used as velocity unit. 
Equation (1161) is one of the main results of this work and 
generalizes the expression found in Ref. [59| to the multi¬ 
band case. 

Equation m, complemented by the coefficients Cj 
and dj derived following the procedure sketched in 
Appendix E allows to obtain the current-phase rela¬ 
tion of the junction. In deriving the scattering coef¬ 
ficients cij, bj, Cj, dj ( j € {a,/3}), Eqs. ©-© have 
to be complemented by the normalization condition, 
JYdx'ijj B (x)^ip B (x) = 1, of the bound state wave 

function ip B (x): T,j= a ,p /-ocJl u j ( x )l 2 + \ v j( x )\ 2 ] dx + 

Ej=a,/3 /o f °°[l«f(*)l 2 + \ v f( x )\ 2 }dx = 1. The normaliza¬ 
tion condition can be satisfied only considering an imagi¬ 
nary part in the expressions of the quasiparticle momenta 
pi and p J h . This accounts for the localized nature of the 
bound state whose wave function cont ains decaying e x- 
ponentials of the type ~ exp[±xr a k F \/l — (E B /A a ) 2 x\ 
and ~ exp[±x r pkF-\/7 2 — (E B /A a ) 2 x\, the decay length 
being comparable with the coherence length £(T) ~ 
Kvf/A a (T). Due to the different decay lengths char¬ 
acterizing the bound state wave function, \E B \ cannot 
exceed A a . Indeed, assuming E B > A a , the quantum 
state is no more localized and its wave function ip B (x) 
is not normalizable. The above arguments explains why 
the bound state energy E B cannot exceed the minimum 
value among the superconducting gaps describing the 
multiband junction. This latter aspect suggests that the 
performances of a multiband Josephson junction can be 
strongly affected by disorder effects and inhomogeneities. 


E b 


2k B T 


M6) 


IV. NUMERICAL RESULTS 


with Uj{x) = 0(—x)Uj(x) + 9(x)uf(x) and analogously 
for Vj(x). Substituting Equation (1141) in the expression 
for J c h, in the absence of spin-sensitive potentials ,5^ ]. 
we obtain (x > 0): 

Jch = Y (x)*d x uf (x)]f(E B ) + 

j 

+ lm[vf(x)d x vf-{x)*](l - f(E B ))}, (15) 

where we explicitly used the thermal equilibrium aver¬ 
ages { 7 t%) = f(E B ) and (%*&) = 1 - f{E B ). Due to 


In this section we provides specific examples of current- 
phase relation derived through Eq. (usd. We assume that 
the temperature dependence of both gaps A^ (j = a, j3) is 

given by A j{T) = Aj (0) tanh ^1.74 *jT c /T — 1^, where 

T c is the critical temperature of the superconducting 
transition. We also assume the validity of the BCS ra¬ 
tio A a (0 )/k B T c = 1.76 between the critical temperature 
and the zero temperature pair potential. Under these as¬ 
sumptions, the parameter 7 = Ap/A a does not depend 
on the temperature. Note that the quantity \o has been 
assigned the value 0.01 in the calculations but that the 
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FIG. 3. Current-phase relations at T = O-OOlTk Panels (a) 
and (c): s± symmetry. Panels (b) and (d): s++ symmetry. 
Parameters are the same of the corresponding panels in Fig.0 
(a) and (b): Z = 0.03, 7 = 3, n = 1.8, r 2 = 1; (c) and (d): 
Z = 0.03, 7 = 3, 77 = 1, r 2 = 1.8. The value of xo has been 
set to xo = 0.01. The various curves are obtained by letting 
parameter s a to vary in the range (0.3, 1) in steps of 0.1, while 
sp = s = 1 . 






T/T c T He 


FIG. 4. Josephson critical current as a function of the re¬ 
duced temperature T/T c . Panels (a) and (c): s± symmetry. 
Panels (b) and (d): s++ symmetry. Parameters are the same 
of the corresponding panels in Fig. [2 ] Z = 0.01, 7 = 3. The 
value of xo has been set to xo = 0.01. (a) and (b) r 1 = 1.8, 
T2 = 1. (c) and (d), r 1 = 1, r2 = 1.8. The critical current 
is in units of eA a (0)/h. The various curves are obtained by 
letting parameter s a to vary in the range (0.3, 1) in steps of 
0.1, while sp = s = 1 . 


results do not depend on the specific value of X 0 ; pro¬ 
vided that xo <C 1 . 

In Figure [3] we present the current-phase relation of two 
junction configurations whose bound states energy has 
been shown in Fig. [2j In particular, in panels (a) and 
(b), we report the current-phase relation computed for 
the s± and the s++ symmetry, respectively, setting the 
model parameters as follows: Z = 0.03, 7 = 3, ry = 1.8, 
7*2 = 1. The different curves are obtained by letting 
the parameter s a to vary in the range (0.3,1) in steps 
of 0.1, while sp = s = 1. One notices a strong sup¬ 
pression of the Josephson current in the s± symmetry 
(panel (a)) compared to the s ++ case (panel (b)). This 
phenomenon originates from the destructive interference 
effects of quasiparticles experiencing two opposite gap 
signs. The destructive interference lowers the Andreev 
reflection probability such that, in perfectly symmetric 
junctions (i.e. r a = rp, s a = sp = s = 1 ), the cur¬ 
rent exactly vanishes. Thus S++ Josephson junctions ex¬ 
hibit greater critical current values (7 C (T —> 0) ~ eA p/K) 
compared to those expected for the s± case (/ C (T —> 0) ~ 
eA a /h). On the other hand, the peculiar functional form 
of the current-phase relation of the s++ symmetry (panel 
(b)) indicates a certain fragility of this junction against 
bias current fluctuations. Indeed, a small bias fluctuation 
(less than few percent of the critical current) can cause, 
at low bias, a relevant phase fluctuation. This stochastic 
phase jump may represent a relevant source of voltage 
fluctuation across the junction which can also drive the 
system towards an ohmic regime. 

Figures (3l(c) and(3](d) are obtained by setting r a = 1.8 
and rp = 1 , while maintaining the remaining parameters 
at the same values fixed in the upper panels of Figure 
[3] In particular, panel (c) shows an s± current-phase 
relation that, as also evident from the Andreev levels 
(Figure [2] (c)), undergoes a 7 r-shift as the s a parameter 
varies between the values 0.3 and 1. The 7 r-shift is absent 
in panel (d) where the current-phase curve is shown for 

the S-|_|_ symmetry. A 7 r-junction presents a free energy 

F(ip) = $o/27r Jq d6I(9) minimum at <p = tt rather than 
at ip = 0. The existence of 7 r-shifted junctions can be ex¬ 
perimentally proved employing SQUID’s measurements. 
While the current-phase relation cannot be simply mea¬ 
sured, the measurement of the temperature dependence 
of the junction critical current is relatively easier. A com¬ 
parison with the experiments of this dependence repre¬ 
sents very often an important tool contributing to identi¬ 
fying the nature of the superconducting pairing. Thus, in 
order to compare our results with the experimental find¬ 
ings, we have derived the temperature dependence of the 
maximum Josephson current (critical current). Figure [4] 
shows the temperature dependence of the critical current 
calculated for the s±- ((a) and (c)) and s ++ -symmetry 
((b) and (d)) by setting the model parameters as done in 
Figs. [2]and[3l Negative critical current values, like those 
presented in panel (c), indicates a 7 r-shifted Josephson 
junction. A common feature of all the panels in Figure 
U is the violation of the Ambegaokar-Baratoff relation, 
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which is usually reported in multiband S/I/S Josephson 
junctions (60). The entity of such violation depends on 
the junction parameters and is more pronounced for the 
s± symmetry (see panels (a) and (c)), where the critical 
current vs temperature curves develop a peculiar posi¬ 
tive curvature. Panels (b) and (d) are described by a de¬ 
formed Ambegaokar-Baratoff relation and suggest that 
a 7 r junction cannot be observed for the s ++ symmetry. 
In this respect we observe that the 0 — 7r transition re¬ 
ported in the present work (see Figure 0(c) and [If c)) is 
the result of a different mechanism compared to the one 
described in Ref. [25[. There the transition originates 
from the competition of two Andreev bound states car¬ 
rying opposite current. As the temperature is lowered 
the state with higher energy eigenvalue is progressively 
excluded from the transport and a sign reversal of the 
critical current is observed. Here we predict the presence 
of a single electron-like bound state and thus a temper¬ 
ature activated 0 — 7r transition cannot be observed as 
an intrinsic effect. However, direct tunneling (a —> a or 
/3 —> f3) and crossed tunneling (a —> ft or /3 —> a) effects 
provide supercurrent contribution with opposite sign and 
thus, depending on the relative strength of these contri¬ 
butions, a 7r-junction can be formed only assuming the s± 
symmetry. The relative strength of the direct and crossed 
tunneling contributions is controlled in the model by the 
overlap parameters s a , sp, and s which are characteristic 
of the interface. According to these arguments, a 0 — 7r 
transition can occur by varying at least one overlap pa¬ 
rameter. This can be done in different ways: (i) changing 
the temperature can induce a lattice deformation at the 
interface which is relevant in determining the overlap pa¬ 
rameters; (ii) experiments performed under controllable 
pressure or strain allow to control lattice distortion and 
thus the overlap parameters. Both the methods (i) and 
(ii) can be used to induce the 0 — 7r transition in multi¬ 
band Josephson junctions. The mechanism of formation 
of a 7r-junction described above can be also recovered in 
the framework of a semiclassical theory of the s±/I/s± 
Josephson junction as developed in Ref. [61]. 

However, as clearly evidenced in Figs. Ela)-(b), the 
0 — 7r transition is strongly affected by the bands ef¬ 
fective mass and can remain unobserved also for an s± 
junction. Under this condition, the magnetic response 
of the system (see Appendix [B] for details) provides an 
indirect probe of the harmonic content of the current- 
phase relation and can help in discriminating the pairing 
symmetry of the junction. In particular, the zero-field 
critical current J c (0) of the s ++ case (see Fig. 0(b)) is 
greater than the corresponding value obtained for the s± 
case (see Fig. 0(a) ). As a consequence, one expects that 
the harmonic content of the current-phase relation of the 
s++ case strongly affects the low-field magnetic diffrac¬ 
tion pattern of the junction. On the other hand, for 
the s± case, the high harmonic contribution is expected 
to be less important due to the lower value of zero-field 
critical current. The above arguments are confirmed by 
Fig. [5] where the normalized critical current J c (<i>)// c (0) 



FIG. 5. Normalized critical current 7 C ($)// C (0) as a function 
of the applied magnetic flux $/$o computed for the s± (panel 
(a)) and the s++ symmetry (panel (b)). Parameters are the 
same of the corresponding panels in Fig. 3: Z = 0.03, 7 = 
3, ri = 1.8, ri = 1, s a = sp — s = 1, xo = 0.01 and 
T = 0.00127. The dashed curve in each panel represents the 
Fraunhofer pattern for comparison. 


is studied as a function of the applied magnetic flux <h/<E>o 
(normalized to the flux quantum $ 0 ) for the s± (panel 
(a)) and the s++ (panel (b)) symmetry, while the model 
parameters are fixed as done for the corresponding pan¬ 
els in Fig. 3. Interestingly, the s++ Josephson junction 
evidences a critical current halving already for flux val¬ 
ues of ~ 0.4 $ Oj while the critical current of the s± case 
is less sensitive to magnetic field effects. 


V. CONCLUSIONS 


We have formulated a minimal model of the dc Joseph¬ 
son effect for multiband superconductors based on the 
quantum waveguide approach. The method is based 
on the analogy between multiband superconductors and 
multibranch networks recently suggested in Ref. 0. 
Accordingly, the subgap bound states wave functions 
(electron- and hole-like), relevant in describing the quan¬ 
tum transport within the short junction limit, are de¬ 
composed using the eigenstates of local branch Hamil¬ 
tonians and the coefficients of such a decomposition are 
found imposing generalized boundary conditions on the 
wave functions. The boundary conditions, never used be¬ 
fore for this problem, are direct generalization of those 
used in Ref. [20} and incorporate a local scattering po¬ 
tential at the interface (Z parameter) and band overlap 
factors s a , sp, and s which define the weight of each 
band in the quasiparticles transport. This provides an 
effective parametrization of the interface effects allowing 
to describe Josephson junctions ranging from the metallic 
{Z <C 1) to the tunnel limit (Z> 1). The bulk properties 
of the superconducting bands are introduced using differ¬ 
ent effective masses (r 3 parameters) which can be deter¬ 
mined by complementary experiments. We have solved 
the scattering problem and we have determined the An¬ 
dreev bound states spectrum and the normalized eigen¬ 
functions. The Josephson current flowing through the 





junction has been computed using a second quantization 
approach which correctly reproduces results obtained us¬ 
ing the phase derivative of the Andreev bound states 
spectrum formula. The second quantization method used 
in the derivation of the Josephson current generalizes the 
result presented in Ref. |59| for a single band supercon¬ 
ductor to the multiband case and represents one of the 
main results of this work. To provide a specific and rele¬ 
vant application of the theory, we have focused our treat¬ 
ment on FeBS-FeBS tunnel junctions. In particular, we 
have derived the current-phase relations and the critical 
currents of symmetric FeBS junctions modeled as super¬ 
conducting systems with two relevant bands. We have 
analyzed s++ and s± symmetries of the FeBS and we 
have shown that a 7r-junction can be observed only for 
s-t symmetry and under appropriate interface conditions, 
which are carefully discussed in the text. Further peculiar 
aspects of FeBS-FeBS tunnel junctions are also discussed. 
The above findings are relevant for the Josephson effect 


theory in multiband systems and can contribute to the 
current debate about the order parameter symmetry of 
iron-based materials. 
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Appendix A: Determination of the scattering 
coefficients 

The matching conditions Eqs. (|8 j) -(|9 jl provide the lin¬ 
ear homogeneous system M • X = 0, where X = 
(a a ,..., d a ,ap ,..., dpY and M is the following matrix 


( 


V 
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(Al) 


The homogeneous system admits non-trivial solution pro¬ 
vided that the condition det( M) = 0 (Eq. (fTOl) of the 
main text) is fulfilled. Since M is not a maximal rank 
matrix, the normalization condition of the bound state 
wave function provides a further condition to determine 
the coefficients vector X. 


Appendix B: Magnetic field dependence of the 
critical current 

We describe the response of a Josephson junction un¬ 
der the application of a sufficiently weak external mag¬ 
netic field H = Hz parallel to the z-axis. We explic¬ 
itly assume the small junction limit in which the trans¬ 
verse dimension L of the junction (parallel to the y-axis) 
is comparable or smaller than the Josephson penetra¬ 
tion length Xj (L < A ? ) [(35}. The junction region, lo¬ 
cated at x = 0, presents a reduced pairing potential 
and thus experiences the maximum magnetic field value, 
while inside the electrodes the polarizing effect of the 
field is effectively screened by supercurrents. Due to this, 
the spatial dependence of the magnetic field is given by 
H = Hz exp(—|x|/A). Since we are interested in the bulk 


effect of H , observing that H ■ zdx = 2A H, we can 
approximate the spatial dependence of the field accord¬ 
ing to the expression H = zj^5(x), where we introduced 
the magnetic flux $ = 2A LH induced by the external 
field and the Dirac delta function S(x). The magnetic 
field H = V x A can be expressed in terms of the vec¬ 
tor potential A = (— ^-S(x), 0, 0) which is not affected 
by self-fields in the considered limit. Assuming that the 
Zeeman term is effectively screened in the bulk of the 
electrodes, the presence of H affects the BdG branch 
Hamiltonian H^ dG (Eq. [T]) only through the substitu¬ 
tion p x —> p x — eA x . The vector potential A x can be 
gauged away by means of the unitary transformation 


U = exp(i7r-^-y6(x)a z ), (Bl) 


with d z a Pauli matrix acting on the particle-hole space 
and $ 0 = h/(2\e\) the elementary flux quantum. The 
transformed branch Hamiltonian Hg\ G = U^Hg dG U un¬ 
der the action of U can be obtained with the following 
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substitutions: 


^ 2 a 2 

-U* 22 H*{x-A x )U 2 2 -)■ 2 ^ 
Aj(x)U^ 1 U 22 -»■ Aj(x), 


(B2) 


~(j) ~ 

while the BdG state associated to H)y d c is given by T, = 
t/ 'T,-. The pair potential Aj(x) takes the same math¬ 
ematical structure of the one presented in Eq. ©, the 
only action of the unitary transformation being the sub¬ 
stitution of the phase difference ip with <p = ip — 

As a consequence the phase difference between the two 
sides of the junction is modulated along the junction and 
this modulation is the source of the magnetic diffraction 
pattern affecting the critical current of the junction. Un¬ 
der our assumptions, the y-dependence of the pair po¬ 
tential is adiabatic compared to the microscopic scale of 
the problem. This statement can be rigorously proved 
using the two-scale perturbation theory [62]. In partic¬ 
ular the wave vector, q = modulating the super¬ 

conducting phase is much smailcr than the Fermi wave 
vectors and thus the phase modulation along the 

y-direction enters only parametrically in the one dimen¬ 
sional problem described in the main text. The validity 
of these arguments also requires that the Cooper pairs 
tunneling with normal incidence represents the dominant 
microscopic process and thus the energy associated with 
transverse modes is negligible compared to the Fermi en¬ 
ergy. When the energy of the transverse modes becomes 
relevant (e.g. for very short junctions with f <C A) on 
the Fermi energy scale, a full treatment of the transverse 
degrees of freedom is required. Hereafter, we focus our 
attention on the case of adiabatic phase variation of <p 
along the transverse dimension of the junction. Under 
this assumption, once the current-phase relation Ij((p) 
has been obtained according to the procedure given in 
the main text, the magnetic field dependence of the crit¬ 
ical current of the junction is obtained according to the 
formula: 


j.L/2 

/ c (<f>) = max / 

¥>e[0,2 tt] J_ L / 2 


dy 


Ij [<p- 2n 


_5_ y_ 
$0 L 


(B3) 


where the current-phase relation presents a parametric 
dependence on y. A further progress can be done ob¬ 
serving that the current-phase relation is an odd func¬ 
tion of the phase difference and thus it can be written as 
Ij(tp) = SriLo a nSm(mp). As a consequence, Eq. (libit) 
takes the form: 

“ sin(n^) 

/ c ($) = max sin (n<p) -(B4) 

■^[o,2U n 


where the coefficients a n can be directly extracted by 
using the relation: 

1 f 27T 

a n = — Ij(ip)sm(rup)dip, (B5) 

n Jo 


where we explicitly used the orthogonality condition 
f 0 n sin(nyj) sin (m<p)d<p = 7rS ntm . Equation (IB4I) provides 
the Fraunhofer diffraction pattern, 

( 7t<£> \ 

Ic(*)/Ic( 0) = 

when a single-harmonic current-phase dependence 
Ij(tp) = ai sin(t^) is considered, while deviations are ex¬ 
pected if the high-harmonic contribution is not negligible. 
For the above reasons, the magnetic diffraction pattern 
J C (<I>) is an indirect probe of the harmonic content of the 
current-phase relation of the junction. Our treatment of 
the magnetic field effects on the junction justifies the va¬ 
lidity of the approach proposed in Refs. HU and the 
classical argument given in Ref. Itul . 
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